Transcranial direct current stimulation of the right frontal eye field in a prosaccade task

R notebook for statistical analysis of the sacc-tDCS dataset. Previous processing:

• Raw data were parsed into events (saccades, fixations, etc.) by the EyeLink data were collected on.
• Events were extracted and saccade measures were computed with a MATLAB script.

# Load some libraries
library("dplyr") # wrangling data frame data
library("tidyr") # tidying data frames
library("ggplot2") # plotting
library("ez") # ANOVA
library("knitr") # R markdown output (html, pdf, etc.)
library("BayesFactor") # Bayesian statistics
library("retimes") # Ex-Gaussian reaction time distributions

Questionnaires

PANAS

Load data

In each session of the experiment, participants completed the Positive and Negative Affect Scale (PANAS) twice:

1. Pre-measurement: Before starting setup of the tDCS and eye tracker
2. Post-measurement: After the electrodes were removed following task completion

As most participants were native Dutch speakers, we also used a Dutch language version of the PANAS, as reported by Engelen et al. (2006).

# Load the data frame
dataFile <- file.path("data", "PANAS.csv")
panasData <- read.csv(dataFile, header = TRUE, sep = ";")
panasData$time <- factor(panasData$time, levels = c("pre","post")) # reorder levels chronologically instead of alphabetically
head(panasData) # show data frame

Factors:

• subject: subject ID (S01, S02, etc)
• session: Whether data are from the first or second session
• stimulation: Whether data are from the anodal or cathodal session
• time: Whether data are from the pre or post measurement

DATA:

The PANAS contains 20 items, split equally into positive and negative affect. Each item is rated on a Likert scale:

1. very slightly or not at all
2. a little
3. moderately
4. quite a bit
5. extremely

To obtain a positive affect score, we sum items: 1 (interested), 3 (excited), 5 (strong), 9 (enthusiastic), 10 (proud), 12 (alert), 14 (inspired), 16 (determined), 17 (attentive) and 19 (active).

To obtain a negative affect score, we sum items: 2 (distressed), 4 (upset), 6 (guilty), 7 (scared), 8 (hostile), 11 (irritable), 13 (ashamed), 15 (nervous), 18 (jittery) and 20 (afraid).

Post-pre differences per item

panasPrePost <- panasData %>%
  group_by(subject,session,stimulation) %>% # for each subject, session, stimulation combination
  select(-time) %>% # don't do subtraction for this column
  summarise_each(funs(diff(.))) # subtract pre and post

Let’s plot the post-pre differences for each item, split for stimulation session.

panasPrePost %>%
  gather(item, score, pos.1.interested:neg.20.afraid) %>% # gather all PANAS columns so it can be used as a factor
  ggplot(aes(item, score)) +
      stat_summary(fun.y = mean, geom = "point", position = position_dodge(width = 0.6), aes(color = stimulation)) +
      geom_hline(yintercept = 0, linetype = "dashed") +
      coord_cartesian(ylim = c(-1,1)) +
      theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

Most scores become negative, meaning the post score was smaller than pre. So they generally feel “less” of everything. This seems particularly the case for positive affect, although participants do feel slightly more proud, determined, and strong.

They become particularly less jittery, nervous, irritable, and also less attentive, active, and excited.

The biggest differences between the sessions are in the positive items, except for irritable.

For some positive items, anodal stimulation results in a more positive score, while cathodal reesults in a more negative score: alert, inspired, determined, and strong.

Composite positive/negative scores

panasComposite <- panasData %>% 
  mutate(positive = rowSums(select(., contains("pos")))) %>% # sum all the positive scores together
  mutate(negative = rowSums(select(., contains("neg")))) %>% # sum all the negative scores together
  select(subject:time, positive, negative) %>% # drop the original PANAS columns
  gather(affect, score, positive, negative) # gather affect so it can be used as a factor

ggplot(panasComposite, aes(time,score, color = stimulation)) +
facet_wrap(~affect) +
  geom_point(position = position_jitter(width = 0.5)) +
  stat_summary(fun.data = mean_cl_normal, position = position_dodge(width = 1))

The negative scores are much lower than the positive overall.

The pre to post differences are larger for positive scores, particularly for the cathodal session.

Statistics: positive scores

Repeated measures ANOVA with factors STIMULATION (anodal vs. cathodal) and TIME (pre vs. post)

panasCompositePositive <- filter(panasComposite, affect == "positive") # new data frame with only positive scores
modelPositive <- ezANOVA(data = data.frame(panasCompositePositive), # Repeated over subjects; type 3 sums of squares (cf. SPSS)
                        dv = .(score), wid = .(subject), within = .(stimulation, time), type = 3)
kable(modelPositive)

	Effect	DFn	DFd	F	p	p<.05	ges
2	stimulation	1	14	0.9180942	0.3542315		0.0080119
3	time	1	14	10.2831973	0.0063350	*	0.0496325
4	stimulation:time	1	14	7.3652913	0.0167928	*	0.0168597

Subjects’ mood becomes less positive after stimulation, particularly in the cathodal session

Statistics: negative scores

Repeated measures ANOVA with factors STIMULATION (anodal vs. cathodal) and TIME (pre vs. post)

panasCompositeNegative <- filter(panasComposite, affect == "negative") # new data frame with only negative scores
modelNegative <- ezANOVA(data = data.frame(panasCompositeNegative), # Repeated over subjects; type 3 sums of squares (cf. SPSS)
                        dv = .(score), wid = .(subject), within = .(stimulation, time), type = 3)
kable(modelNegative)

	Effect	DFn	DFd	F	p	p<.05	ges
2	stimulation	1	14	0.1054217	0.7502187		0.0007105
3	time	1	14	16.2682403	0.0012327	*	0.0845895
4	stimulation:time	1	14	3.8757669	0.0691122		0.0101630

Subjects’ mood also become less negative after stimulation…

tDCS sensations

After each session, subjects also completed a questionnaire probing tDCS sensations.

Load data

# Load the data frame
dataFile <- file.path("data", "tdcs_sensations.csv")
sensData <- read.csv(dataFile, header = TRUE, sep = ";", na.strings = "")
head(sensData) # show data frame

They were asked to which degree the following sensations were present during stimulation: tingling, itching sensation, burning sensation, pain, headache, fatigue, dizziness and nausea. Each was rated on a scale from 0-4:

0. none
1. a little
2. moderate
3. strong
4. very strong

They also rated their confidence that the sensations were caused by the stimulation on a scale from 0-4 (columns starting with conf.):

0. n/a (meaning they rated the sensation a 0 on the previous scale)
1. unlikely
2. possibly
3. likely
4. very likely

Finally, they filled in whether they felt one electrode more than the other (felt.more); (and if so, which and for which sensations).

Factors:

• subject: subject ID (S01, S02, etc)
• session: Whether data are from the first or second session
• stimulation: Whether data are from the anodal or cathodal session

Plot sensation distribution

sensData %>%
  select(everything(), -contains("conf"), -felt.more) %>% # drop the confidence columns
  gather(sensation, rating, itching:nausea) %>% # gather sensations so they can be used as a factor
  ggplot(aes(rating, fill = stimulation)) +
    facet_wrap(~sensation) +
    geom_histogram(position = "dodge", binwidth = 0.5) +
    xlim(0.5,4.5) +
    ggtitle("Sensations over 15 anodal sessions, 15 cathodal sessions")

Nausea was never experienced; dizziness and pain are also rare, and mild.

Burning, itching and tingling are most frequently experienced and to a higher degree.

Plot confidence distributions

sensData %>%
  select(contains("conf"), subject, session, stimulation) %>% # drop the rating columns
  gather(sensation, rating, conf.itching:conf.nausea) %>%
  ggplot(aes(rating, fill = stimulation)) +
    facet_wrap(~sensation) +
    geom_histogram(position = "dodge", binwidth = 0.5) +
    xlim(0.5,4.5) +
    ggtitle("Confidence over 15 anodal sessions, 15 cathodal sessions")

For “local” sensations like burning, tingling and dizziness, subjects have high confidence that these are due to tDCS. For more diffuse sensations, like fatigue and headache, ratings are very low, so their occurence might just be due to performing the task for an extended period of time.

Sensation difference between anode and cathode

sensData %>%
  select(felt.more, subject, session, stimulation) %>% # keep only the relevant column
  mutate(felt.more = factor(felt.more, levels = c("anode", "equal", "cathode"))) %>% # keep colors consistent
  ggplot(aes(stimulation, fill = felt.more)) +
    geom_bar() +
  scale_y_continuous(breaks = 1:15)

There’s a pretty even split between which electrode people feel the most, so there don’t appear to be worrisome biases.

During anodal stimulation, the anode is felt more than the cathode; during cathodal stimulation the cathode is felt more than the anode. In other words, the electrode over the FEF is always felt more than the forehead electrode.

About a third of people did not indicate they felt one more than the other; this is slightly more in the cathodal session.

Statistics

We will test for every sensation separately whether there was a difference between the anodal and cathodal sessions. Mann-Whitney-U tests are most appropriate here, as the data are ordinal (Likert) and do not look normally distributed.

Itching

wilcox.test(sensData$itching[sensData$stimulation == "anodal"], sensData$itching[sensData$stimulation == "cathodal"])

cannot compute exact p-value with ties

    Wilcoxon rank sum test with continuity correction

data:  sensData$itching[sensData$stimulation == "anodal"] and sensData$itching[sensData$stimulation == "cathodal"]
W = 122, p-value = 0.6828
alternative hypothesis: true location shift is not equal to 0

Tingling

wilcox.test(sensData$tingling[sensData$stimulation == "anodal"], sensData$tingling[sensData$stimulation == "cathodal"])

cannot compute exact p-value with ties

    Wilcoxon rank sum test with continuity correction

data:  sensData$tingling[sensData$stimulation == "anodal"] and sensData$tingling[sensData$stimulation == "cathodal"]
W = 93, p-value = 0.4097
alternative hypothesis: true location shift is not equal to 0

Burning

wilcox.test(sensData$burning[sensData$stimulation == "anodal"], sensData$burning[sensData$stimulation == "cathodal"])

cannot compute exact p-value with ties

    Wilcoxon rank sum test with continuity correction

data:  sensData$burning[sensData$stimulation == "anodal"] and sensData$burning[sensData$stimulation == "cathodal"]
W = 79.5, p-value = 0.1571
alternative hypothesis: true location shift is not equal to 0

Pain

wilcox.test(sensData$pain[sensData$stimulation == "anodal"], sensData$pain[sensData$stimulation == "cathodal"])

cannot compute exact p-value with ties

    Wilcoxon rank sum test with continuity correction

data:  sensData$pain[sensData$stimulation == "anodal"] and sensData$pain[sensData$stimulation == "cathodal"]
W = 94, p-value = 0.3126
alternative hypothesis: true location shift is not equal to 0

Headache

wilcox.test(sensData$headache[sensData$stimulation == "anodal"], sensData$headache[sensData$stimulation == "cathodal"])

cannot compute exact p-value with ties

    Wilcoxon rank sum test with continuity correction

data:  sensData$headache[sensData$stimulation == "anodal"] and sensData$headache[sensData$stimulation == "cathodal"]
W = 118, p-value = 0.7782
alternative hypothesis: true location shift is not equal to 0

Fatigue

wilcox.test(sensData$fatigue[sensData$stimulation == "anodal"], sensData$fatigue[sensData$stimulation == "cathodal"])

cannot compute exact p-value with ties

    Wilcoxon rank sum test with continuity correction

data:  sensData$fatigue[sensData$stimulation == "anodal"] and sensData$fatigue[sensData$stimulation == "cathodal"]
W = 100, p-value = 0.5757
alternative hypothesis: true location shift is not equal to 0

Dizziness

wilcox.test(sensData$dizziness[sensData$stimulation == "anodal"], sensData$dizziness[sensData$stimulation == "cathodal"])

cannot compute exact p-value with ties

    Wilcoxon rank sum test with continuity correction

data:  sensData$dizziness[sensData$stimulation == "anodal"] and sensData$dizziness[sensData$stimulation == "cathodal"]
W = 105.5, p-value = 0.6045
alternative hypothesis: true location shift is not equal to 0

Nausea

wilcox.test(sensData$nausea[sensData$stimulation == "anodal"], sensData$nausea[sensData$stimulation == "cathodal"])

cannot compute exact p-value with ties

    Wilcoxon rank sum test with continuity correction

data:  sensData$nausea[sensData$stimulation == "anodal"] and sensData$nausea[sensData$stimulation == "cathodal"]
W = 112.5, p-value = NA
alternative hypothesis: true location shift is not equal to 0

Nausea was rated 0 by everyone for both anodal and cathodal.

Load eye data

The .csv file with the eye tracking data was created in MATLAB.

# Load the data frame
dataFile <- file.path("data", "sacc-tDCS_data.csv")
groupData <- read.csv(dataFile, header = TRUE, na.strings = "NaN")

# Format data frame
groupData <- select(groupData, -X) # drop final empty column named "X"
groupData$leg <- factor(groupData$leg, levels = c("pre", "tDCS", "post"))

kable(head(groupData))

subject	stimulation	leg	block	trial	type	direction	deviation.start	deviation.end	amplitude	latency
S01	anodal	pre	1	1	lateral	right	0.462897	0.170648	8.02463	433
S01	anodal	pre	1	1	center	left	0.459092	1.035020	7.16262	439
S01	anodal	pre	1	2	lateral	right	0.344561	0.409786	7.74873	291
S01	anodal	pre	1	2	center	left	0.550230	0.505107	7.43233	198
S01	anodal	pre	1	3	lateral	right	0.514736	0.612098	7.61080	281
S01	anodal	pre	1	3	center	left	0.620728	1.627990	6.35043	376

• subject: subject IDj
• stimulation: Whether data are from the anodal or cathodal session
• leg: Whether data are before (pre), during (tDCS), or after (post) tDCS
• block: After each block participant had a brief break and tracker was recalibrated
• trial: trial number within a block
• type:
  • lateral - fixation in center of display, saccade made towards the periphery
  • center - fixation in periphery, saccade made back towards the center of the display
• direction: left for saccades towards the left of current fixation position; right for saccades to the right
• deviation.start : distance (in visual angle) from saccade start point to fixation
• deviation.end: distance (in visual angle) from saccade end point to target location
• amplitude: distance (in visual angle) between saccade start and end point
• latency: time (in ms) from target onset to start of saccade

Basic inspection

Reaction time distributions

Histograms for each subject

histType <- ggplot(groupData, aes(latency, fill = type)) +
  facet_wrap(~subject, nrow = 3, scales ="free_y") +
  geom_histogram(binwidth = 5, color = "grey50", size = .2) +
  xlim(-50,300)
histType

Stray observations:

• Center saccades are much faster than lateral, though not to the same degree in all subjects
• Some have a fat short latency tail - too fast saccades that are virtually all towards the center: S05, S10, S11
• Some appear bimodal, but when split for type this is generally because center saccades are faster: S05, S06, S11
• Some are super sharp (S8); others really broad (S01)
• S01 has a very typical looking distribution, but slow

Stimulation effects across subjects

dens <- ggplot(groupData, aes(latency, color = stimulation, linetype = leg)) +
  facet_grid(type ~ direction) +
  geom_density() +
  xlim(0, 250) +
  scale_color_brewer(palette = "Set1")
dens

Anodal vs. cathodal in each subject

denstDCS <- ggplot(groupData[groupData$leg == 'tDCS' & groupData$type == "lateral", ], aes(latency, color = stimulation)) +
  facet_wrap(~subject, nrow = 3, scales ="free_y") +
  geom_density() +
  xlim(0, 250) +
  scale_color_brewer(palette = "Set1") +
  ggtitle('Lateral saccades, tDCS block')
denstDCS

Anodal session in each subject

densLegAnodal <- ggplot(groupData[groupData$stimulation == 'anodal' & groupData$type == "lateral", ], aes(latency, color = leg)) +
  facet_wrap(~subject, nrow = 3, scales ="free_y") +
  geom_density() +
  xlim(0, 250) +
  ggtitle('Lateral saccades, anodal')
densLegAnodal

Cathodal session in each subject

densLegCathodal <- ggplot(groupData[groupData$stimulation == 'cathodal' & groupData$type == "lateral", ], aes(latency, color = leg)) +
  facet_wrap(~subject, nrow = 3, scales ="free_y") +
  geom_density() +
  xlim(0, 250) +
  ggtitle('Lateral saccades, cathodal')
densLegCathodal

Outliers

Mark outliers

Criteria for outliers:

• Discard fast saccades, with a latency of 50 ms or less
• Discard slow saccades, saccades with a latency of 400 ms or more
• Discard inaccurate fixations, with saccade starting point more than 1.8 degrees or more away from fixation
• Discard faulty saccades, with saccade end point more than 8 degree or more away from the target

In Kanai et al. (2012), this was:

• Fast saccades: 50 ms
• Slow saccades: 400 ms
• Bad fixations: 1.8 degrees
• Faulty saccades: opposite hemifield of target (here, that would be 8 degrees as targets were that eccentric)

# Mark outliers
groupData <- mutate(groupData, outlier = FALSE, # fill vector with FALSE for all trials
                    outlier = ifelse(latency < tooFast, "fast", outlier), # mark too fast trials as "fast"
                    outlier = ifelse(latency > tooSlow, "slow", outlier), # mark too slow trials as "slow"
                    outlier = ifelse(deviation.start > badFix, "fixation", outlier), # mark bad fixations as "fixation"
                    outlier = ifelse(deviation.end > badSacc, "saccade", outlier)) # mark inaccurate saccades as "saccade"
groupDataClean <- filter(groupData, outlier == FALSE) # make new data frame without outlier trials

Plot outliers per subject

outlierPlot <- ggplot(groupData[groupData$outlier != FALSE, ], aes(interaction(stimulation,leg,block,trial), latency, color = outlier, shape = type)) +
  facet_wrap(~subject, nrow = 4, scales = "free_y") +
  geom_point() +
  geom_hline(yintercept = 0) +
  geom_hline(yintercept = tooFast, linetype = "dashed") +
  geom_hline(yintercept = tooSlow, linetype = "dashed") +
  xlab('Trial') +
  theme(axis.text.x=element_blank(), # remove x-axis (just trial count)
  axis.ticks.x=element_blank())
outlierPlot

outlierCount <- groupData %>%
  group_by(subject) %>%
  summarize(outlier_count = sum(outlier != FALSE, na.rm = TRUE))
kable(outlierCount, caption = "Number of outlier trials per subject")

subject	outlier_count
S01	180
S02	35
S03	282
S04	154
S05	536
S06	432
S07	48
S08	180
S09	116
S10	730
S11	303
S12	101
S13	176
S14	426
S15	155
S16	1026
S17	394
S18	28
S19	289
S20	383

Stray observations:

Differences between subjects:

• Some subjects have barely any outliers (S02, S07)
• Most subjects have quite a few outliers, especially S05,S06, S10 and S1. The mean is 299 out of 5760 trials in total
• Only S01 has a sizable amount of slow saccades
• S14 has a lot of negative latencies (and also long positive latencies)
• Those subjects with many fast saccades also tend to have many bad fixations (S06, S10, but see S08)

General patterns:

• Occurence of outliers seems stable throughout the session: there aren’t more/less in the beginning / end
• There are very few inaccurate saccades (makes sense, because task is easy and criterion is not strict)
• Most outliers are too fast saccades and bad fixations
• Slow saccades are lateral, fast saccades are to the center, because only the latter are predictable
• Bad fixations appear to be mostly center saccades (but this varies a lot). Perhaps the eyes already drift back towards the center, before executing the saccade?
• Or is the source of bad fixations simply poor quality of the eye tracker data? e.g. people are actually fixating, but due to drift it appears they are not.

Median reaction time

Here we simply extract median RTs for each condition and use a repeated measures ANOVA for statistical analysis, following Kanai et al. (2012).

Data per block

# Compute median in each condition
latencyMedian <- groupDataClean %>%
  group_by(subject,stimulation,leg,block,type,direction) %>%
  summarise(latency = median(latency, na.rm = TRUE))

Full factorial plot

# Plot out all the data
fullPlot <- ggplot(latencyMedian, aes(interaction(block,leg), latency, color = stimulation, shape = stimulation)) +
  facet_grid(type ~ direction) +
  stat_summary(fun.y = mean, geom = "point", size = 3) +
  stat_summary(fun.y = mean, geom = "line", aes(group = stimulation), size = 1)
fullPlot

There is a pretty sizable difference between the anodal & cathodal sessions in the baseline already: cathodal is always faster than anodal.

The above plot also shows a lot of variability, so it might be best to average over each 3 consecutive blocks so the data come in 15-minute intervals.

15-minute intervals (collapse 3 blocks)

# Compute median per leg
latencyMedianLeg <- groupDataClean %>%
  group_by(subject,stimulation,direction,type) %>% 
  summarise(baseline = median(latency[leg == "pre"], na.rm = TRUE), # take average of 3 blocks, make new column
            tDCS = median(latency[leg == "tDCS"], na.rm = TRUE),
            post.1 = median(latency[leg == "post" & block <= 3], na.rm = TRUE),
            post.2 = median(latency[leg == "post" & block >= 4], na.rm = TRUE)) %>%
 gather(leg,latency,baseline,tDCS,post.1,post.2) %>% # gather new columns to use as factor
 mutate(leg = factor(leg, levels = c("baseline", "tDCS", "post.1", "post.2"))) # reorder factor levels

Line plot per leg, individual subjects

Now make the same plot, but for the data collapsed over 3 blocks. Also draw the plot for each individual subject, so we can see which subjects drive the baseline difference, and in which direction the stimulation effect goes for each subject (if there is any).

kanaiSubsPlot <- ggplot(latencyMedianLeg, aes(leg, latency, color = stimulation, shape = type)) +         
  facet_wrap(~subject, nrow = 5) +
  stat_summary(fun.y = mean, geom = "point", size = 2) +
  stat_summary(fun.y = mean, geom = "line", aes(group = stimulation))
kanaiSubsPlot

There are a couple of subjects that show large differences in the baseline already, especially S01.

Compute magnitude of baseline difference

Let’s look at the size of the baseline difference per subject.

baselineDiff <- latencyMedianLeg %>% 
  filter(leg == "baseline") %>% # keep only baseline data
  spread(stimulation, latency) %>% # make separate columns for anodal and cathodal
  mutate(latency.diff = anodal - cathodal) %>% # subtract the difference
  group_by(subject) %>% 
  summarise(latency.diff = mean(latency.diff))# keep the average difference per subject
kable(baselineDiff, caption = 'Difference between baseline saccade latencies in anodal and cathodal session')

subject	latency.diff
S01	39.125
S02	18.625
S03	2.375
S04	-6.125
S05	-17.500
S06	-3.875
S07	5.625
S08	11.125
S09	23.750
S10	14.750
S11	1.000
S12	8.375
S13	4.000
S14	-10.375
S15	1.875
S16	-18.750
S17	-5.500
S18	-2.875
S19	1.000
S20	4.750

There are a few subjects with substantial latency differences between sessions; the three largest of which are all slower in the anodal session (positive values). The mean difference is 3.6 ms.

Exclude S01 and S16

It seems reasonable to exclude subject 1 from further analysis, because of three things:

• S01 shows the largest baseline difference between the anodal and cathodal sessions.
• Overall, S01 is much slower than all other subjects, by several tens of milliseconds.
• The first session for S01 was aborted because S01 was not feeling well. It was repeated at a later date, so S01 came to the lab 3 times in total instead of two.

However, S01 also seems to show the largest effect of stimulation, for both anodal and cathodal stimulation. It is even in the expected direction: faster saccades than baseline in the anodal session, slower in cathodal.

S16’s data quality was poor in the first session and also shows a large baseline difference

latencyMedianLegExcl <- filter(latencyMedianLeg, subject != "S01") # discard rows from S01
latencyMedianLegExcl$subject <- factor(latencyMedianLegExcl$subject) #remake factor as it now has one level less

Line plot per leg over all subjects

Now that we’ve excluded subject(s) and collapsed data over blocks, let’s look at the group average plot for the first time.

kanaiPlot <- ggplot(latencyMedianLegExcl, aes(leg, latency, color = stimulation, shape = stimulation)) +         
  facet_grid(type ~ direction) +
  stat_summary(fun.y = mean, geom = "point", size = 3) +
  stat_summary(fun.y = mean, geom = "line", aes(group = stimulation), size = 1) +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", width = 0.3)
kanaiPlot

Even without S01, there is still a small baseline difference of 1.7 ms. There also seems to be a tiny effect in the center-left saccade condition that’s in the expected direction for both stimulation sessions, but it’s about the size of the baseline difference.

Subtract baseline

The baseline difference is not informative and could obscure real changes from baseline within subjects. Therefore, subtract the baseline from each subsequent measurement.

latencyMedianBaseline <- latencyMedianLegExcl %>%
  group_by(subject,stimulation,direction,type) %>% # for each condition, subtract baseline scores and make new columns
  summarise(tDCS = latency[leg == "tDCS"] - latency[leg == "baseline"], 
           post.1 = latency[leg == "post.1"] - latency[leg == "baseline"],
           post.2 = latency[leg == "post.2"] - latency[leg == "baseline"]) %>%
  gather(leg, latency, tDCS, post.1, post.2)  %>% # gather new columns to use as factor 
  mutate(leg = factor(leg, levels = c("tDCS", "post.1", "post.2"))) # reorder factor levels

Line plots per leg from baseline

kanaiPlotBase <- ggplot(latencyMedianBaseline, aes(leg, latency, color = stimulation, shape = stimulation)) +         
  facet_grid(type ~ direction) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  stat_summary(fun.y = mean, geom = "point", size = 3) +
  stat_summary(fun.y = mean, geom = "line", aes(group = stimulation), size = 1) +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", width = 0.3) +
  coord_cartesian(ylim = c(-15,15))
kanaiPlotBase

All changes from baseline are really tiny: 5 ms or less.

Individual subjects, anodal session

kanaiPlotBaseSubsAnodal <- ggplot(latencyMedianBaseline[latencyMedianBaseline$stimulation == "anodal", ], aes(leg, latency)) +
  facet_grid(type ~ direction) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_line(aes(group=subject,color=subject)) +
  stat_summary(fun.y = mean, aes(group = stimulation), geom = "line") +
  stat_summary(fun.y = mean, geom = "point") +
  ggtitle("Anodal difference from baseline")
kanaiPlotBaseSubsAnodal

For the center saccades, most subjects follow the expected direction and are below the line. The left-lateral effect that is of primary interest is erratic though: only ~10 subjects are below the line, and none show an online-effect larger than 12 ms.

Individual subjects, cathodal session

kanaiPlotBaseSubsCathodal <- ggplot(latencyMedianBaseline[latencyMedianBaseline$stimulation == "cathodal", ], aes(leg, latency)) +
  facet_grid(type ~ direction) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_line(aes(group=subject,color=subject)) +
  stat_summary(fun.y = mean, aes(group = stimulation), geom = "line") +
  stat_summary(fun.y = mean, geom = "point") + 
  ggtitle("Cathodal difference from baseline")
kanaiPlotBaseSubsCathodal

All of this is split pretty much 50-50, hence the average difference hovering around 0.

Statistics

T-tests for baseline differences

First, test for baseline differences for each combination of DIRECTION (left, right) and TYPE (center, lateral).

# keep only the baseline for left-lateral-saccades
latencyMedianPreLeftLateral <- latencyMedianLegExcl %>%
  filter(leg == "baseline", type == "lateral", direction == "left") %>%
  select(-leg,-type,-direction) %>%
  spread(stimulation, latency)

Adding missing grouping variables: `direction`

# keep only the baseline for right-lateral-saccades
latencyMedianPreRightLateral <- latencyMedianLegExcl %>%
  filter(leg == "baseline", type == "lateral", direction == "right") %>%
  select(-leg,-type,-direction) %>%
  spread(stimulation, latency)

Adding missing grouping variables: `direction`

# keep only the baseline for right-center-saccades
latencyMedianPreRightCenter <- latencyMedianLegExcl %>%
  filter(leg == "baseline", type == "center", direction == "right") %>%
  select(-leg,-type,-direction) %>%
  spread(stimulation, latency)

Adding missing grouping variables: `direction`

# keep only the baseline for left-center-saccades
latencyMedianPreLeftCenter <- latencyMedianLegExcl %>%
  filter(leg == "baseline", type == "center", direction == "left") %>%
  select(-leg,-type,-direction) %>%
  spread(stimulation, latency)

Adding missing grouping variables: `direction`

t.test(latencyMedianPreLeftLateral$anodal,latencyMedianPreLeftLateral$cathodal, paired = TRUE )

    Paired t-test

data:  latencyMedianPreLeftLateral$anodal and latencyMedianPreLeftLateral$cathodal
t = 0.59438, df = 18, p-value = 0.5597
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -3.668569  6.563306
sample estimates:
mean of the differences 
               1.447368 

t.test(latencyMedianPreRightLateral$anodal,latencyMedianPreRightLateral$cathodal, paired = TRUE )

    Paired t-test

data:  latencyMedianPreRightLateral$anodal and latencyMedianPreRightLateral$cathodal
t = 1.299, df = 18, p-value = 0.2104
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -2.290841  9.711893
sample estimates:
mean of the differences 
               3.710526 

t.test(latencyMedianPreLeftCenter$anodal,latencyMedianPreLeftCenter$cathodal, paired = TRUE )

    Paired t-test

data:  latencyMedianPreLeftCenter$anodal and latencyMedianPreLeftCenter$cathodal
t = 0.22858, df = 18, p-value = 0.8218
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -5.820180  7.241232
sample estimates:
mean of the differences 
              0.7105263 

t.test(latencyMedianPreRightCenter$anodal,latencyMedianPreRightCenter$cathodal, paired = TRUE )

    Paired t-test

data:  latencyMedianPreRightCenter$anodal and latencyMedianPreRightCenter$cathodal
t = 0.28503, df = 18, p-value = 0.7789
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -5.867854  7.709959
sample estimates:
mean of the differences 
              0.9210526 

None of the baseline differences are significant, so there is reason to leave the baseline data in the rest of the analyses.

Omnibus anova - saccade latency

Data: Outliers removed, collapesed into 15-minute intervals.

Dependent measure: saccadic latency

Factors:

• STIMULATION (anodal vs. cathodal)
• LEG (baseline, tDCS, post.1, post.2)
• TYPE (lateral vs. center)
• DIRECTION (left vs. right)

modelOmni <- ezANOVA(data = data.frame(latencyMedianLegExcl), # Repeated over subjects; type 3 sums of squares (cf. SPSS)
                        dv = .(latency), wid = .(subject), within = .(stimulation, leg, type, direction), type = 3)
kable(modelOmni$ANOVA)

	Effect	DFn	DFd	F	p	p<.05	ges
2	stimulation	1	18	0.5448336	0.4699513		0.0026522
3	leg	3	54	0.5719557	0.6359004		0.0010567
4	type	1	18	26.5163705	0.0000673	*	0.2031419
5	direction	1	18	0.0443702	0.8355326		0.0003309
6	stimulation:leg	3	54	0.0276830	0.9937138		0.0000243
7	stimulation:type	1	18	0.7878197	0.3864669		0.0008186
8	leg:type	3	54	2.9395372	0.0412563	*	0.0015450
9	stimulation:direction	1	18	1.0471920	0.3197118		0.0004897
10	leg:direction	3	54	1.1003477	0.3570242		0.0002611
11	type:direction	1	18	1.3802333	0.2553696		0.0013139
12	stimulation:leg:type	3	54	0.0940556	0.9630059		0.0000290
13	stimulation:leg:direction	3	54	4.1363319	0.0103488	*	0.0007042
14	stimulation:type:direction	1	18	1.8347183	0.1923311		0.0001654
15	leg:type:direction	3	54	2.7550809	0.0512305		0.0006159
16	stimulation:leg:type:direction	3	54	0.9293330	0.4329215		0.0001384

Main effect: type

latencyMedianLegExcl %>%
  group_by(subject,type) %>%
  summarise(latency = mean(latency)) %>%
  ggplot(aes(type, latency)) +
  stat_summary(fun.data = mean_cl_normal, size = 1) +
  geom_jitter(width = 0.25)

This simply reflects that center saccades are faster than lateral saccades, because the location of the target is known.

Interaction: leg by type

latencyMedianLegExcl %>%
  group_by(subject,leg,type) %>%
  summarise(latency = mean(latency)) %>%
  ggplot(aes(leg, latency, shape = type)) +
  stat_summary(fun.y = mean, geom = "point") +
  stat_summary(fun.y = mean, geom = "line", aes(group = type, linetype = type))

The effect of TYPE becomes larger over time: center saccades get faster, lateral saccades become slower.

Interaction: leg by type by direction

latencyMedianLegExcl %>%
  group_by(subject,leg,type,direction) %>%
  summarise(latency = mean(latency)) %>%
  ggplot(aes(leg, latency, shape = type)) +
  facet_wrap(~direction) +
  stat_summary(fun.y = mean, geom = "point") +
  stat_summary(fun.y = mean, geom = "line", aes(group = type, linetype = type))

This interaction of LEG and TYPE appears to occur to a lesser extent for left saccades. The “tDCS” block is the deviant here, but there is no interaction with stimulation.

ANOVA matching Kanai et al. (2012) - saccade latency

Differing from the previous omnibus analysis, Kanai et al. (2012) analysed shifts from baseline and only had lateral saccades.

Data:

• Outliers removed
• Collapsed into 15-minute intervals
• Subtract the baseline from each subsequent block
• Discard center, keep only lateral saccades

Dependent measure: saccadic latency

Factors:

• STIMULATION (anodal vs. cathodal)
• LEG (tDCS, post.1, post.2)
• DIRECTION (left vs. right)

latencyMedianBaselineLateral <- filter(latencyMedianBaseline, type == "lateral") # keep only lateral saccades
modelKanai <- ezANOVA(data = data.frame(latencyMedianBaselineLateral),
                        dv = .(latency), wid = .(subject), within = .(stimulation,leg,direction), type = 3)
# OR, without the EZ package:
# modelKanai=aov(latency~stimulation*leg*direction + Error(subject/(stimulation*leg*direction)),data=latencyMedianBaselineLateral)
# summary(modelKanai)
kable(modelKanai$ANOVA)

	Effect	DFn	DFd	F	p	p<.05	ges
2	stimulation	1	18	0.1056396	0.7489113		0.0011603
3	leg	2	36	3.0852605	0.0579807		0.0254120
4	direction	1	18	0.0006761	0.9795419		0.0000049
5	stimulation:leg	2	36	0.0500856	0.9512141		0.0001864
6	stimulation:direction	1	18	0.1124817	0.7412135		0.0002421
7	leg:direction	2	36	1.6749929	0.2015777		0.0014317
8	stimulation:leg:direction	2	36	4.0522338	0.0258677	*	0.0060139

Main effect of direction

latencyMedianBaselineLateral %>%
  group_by(subject,direction) %>%
  summarise(latency = mean(latency)) %>%
  ggplot(aes(direction, latency)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  stat_summary(fun.data = mean_cl_normal, size = 1) +
  geom_jitter(width = 0.25)

Right saccades become slower with respect to the baseline; left saccades do not. However, the effect is tiny (1-2 ms) and there is no interaction with STIMULATION.

ANOVA matching Kanai et al. (2012) - center saccades

Repeat the same ANOVA, but now for center saccades (which Kanai did not have).

latencyMedianBaselineCenter <- filter(latencyMedianBaseline, type == "center") # keep only lateral saccades
modelKanaiCenter <- ezANOVA(data = data.frame(latencyMedianBaselineCenter),
                        dv = .(latency), wid = .(subject), within = .(stimulation,leg,direction), type = 3)
# OR, without the EZ package:
# modelKanai=aov(latency~stimulation*leg*direction + Error(subject/(stimulation*leg*direction)),data=latencyMedianBaselineLateral)
# summary(modelKanai)
kable(modelKanaiCenter$ANOVA)

	Effect	DFn	DFd	F	p	p<.05	ges
2	stimulation	1	18	0.0125981	0.9118744		0.0001639
3	leg	2	36	0.7986263	0.4577571		0.0027491
4	direction	1	18	4.2193901	0.0547812		0.0127626
5	stimulation:leg	2	36	0.0519098	0.9494854		0.0001221
6	stimulation:direction	1	18	0.0980885	0.7577357		0.0001511
7	leg:direction	2	36	1.8879711	0.1660644		0.0025036
8	stimulation:leg:direction	2	36	3.2728364	0.0494365	*	0.0032078

ANOVA without subtracting baseline, for lateral saccades - saccade latency

Kanai et al. analyzed the baseline-subtracted data; now repeat the same ANOVA with the baseline block still present.

Data:

• Outliers removed
• Collapsed into 15-minute intervals
• Discard center, keep only lateral saccades

Dependent measure: saccadic latency

Factors:

• STIMULATION (anodal vs. cathodal)
• LEG (tDCS, post.1, post.2)
• DIRECTION (left vs. right)

# keep only lateral saccades
latencyMedianLateral <- latencyMedianLegExcl %>%
  filter(type == "lateral") %>%
  select(-type)
modelLateral <- ezANOVA(data = data.frame(latencyMedianLateral),
                        dv = .(latency), wid = .(subject), within = .(stimulation, leg, direction), type = 3)
kable(modelLateral$ANOVA)

	Effect	DFn	DFd	F	p	p<.05	ges
2	stimulation	1	18	1.1896030	0.2897954		0.0057760
3	leg	3	54	2.1450309	0.1052250		0.0034252
4	direction	1	18	0.3432586	0.5652290		0.0026728
5	stimulation:leg	3	54	0.0736893	0.9738537		0.0000562
6	stimulation:direction	1	18	1.7519802	0.2021940		0.0011053
7	leg:direction	3	54	0.5348350	0.6603829		0.0001691
8	stimulation:leg:direction	3	54	3.0084971	0.0380553	*	0.0007197

Leg by direction interaction

latencyMedianLateral %>%
  group_by(subject,leg,direction) %>%
  summarise(latency = mean(latency)) %>%
  ggplot(aes(leg, latency, shape = direction)) +
  stat_summary(fun.y = mean, geom = "point") +
  stat_summary(fun.y = mean, geom = "line", aes(group = direction, linetype = direction))

Pretty erratic pattern, and no interaction with STIMULATION either.

ANOVA without subtracting baseline, for center saccades - saccade latency

Same as previous, but now for center saccades.

latencyMedianCenter <- latencyMedianLegExcl %>%
  filter(type == "center") %>%
  select(-type)
modelCenter <- ezANOVA(data = data.frame(latencyMedianCenter),
                        dv = .(latency), wid = .(subject), within = .(
                          stimulation, leg, direction), type = 3)
kable(modelCenter$ANOVA)

	Effect	DFn	DFd	F	p	p<.05	ges
2	stimulation	1	18	0.0818013	0.7781366		0.0005891
3	leg	3	54	0.4720792	0.7030047		0.0015717
4	direction	1	18	0.0399289	0.8438598		0.0003658
5	stimulation:leg	3	54	0.0290943	0.9932360		0.0000499
6	stimulation:direction	1	18	0.2061376	0.6552385		0.0000962
7	leg:direction	3	54	2.7402193	0.0521336		0.0017524
8	stimulation:leg:direction	3	54	2.3794025	0.0797697		0.0009944

Bayesian linear mixed effects matching Kanai - saccade latency

Test against the null model

bfKanaiNull = anovaBF(latency~stimulation*leg*direction+subject, data = data.frame(latencyMedianBaselineLateral), whichModels="withmain", whichRandom = "subject", progress = FALSE, iterations = 100000) # compute Bayes Factors
bfKanaiNull = sort(bfKanaiNull, decreasing = TRUE) # sort such that winning model is at the top

First we compare all models to the most simple (null) model, which is the intercept only + random effect model: latency ~ subject. This does not test for effects of SUBJECT but models it as a nuisance factor (whichRandom = "subject"). In addition, to decrease the model space, we do not consider models that have an interaction without the corresponding main effects (whichModels = "withmain").

kable(select(extractBF(bfKanaiNull), bf)) # show only the Bayes factors in a table

	bf
leg + subject	1.9394726
stimulation + leg + subject	0.3327031
direction + leg + subject	0.2801473
stimulation + subject	0.1716592
direction + subject	0.1427363
stimulation + direction + leg + subject	0.0479931
stimulation + leg + stimulation:leg + subject	0.0279650
direction + leg + direction:leg + subject	0.0273854
stimulation + direction + subject	0.0243691
stimulation + direction + stimulation:direction + leg + subject	0.0095508
stimulation + direction + stimulation:direction + subject	0.0052792
stimulation + direction + leg + direction:leg + subject	0.0047485
stimulation + direction + leg + stimulation:leg + subject	0.0039681
stimulation + direction + stimulation:direction + leg + direction:leg + subject	0.0009712
stimulation + direction + stimulation:direction + leg + stimulation:leg + subject	0.0008092
stimulation + direction + leg + stimulation:leg + direction:leg + subject	0.0003951
stimulation + direction + stimulation:direction + leg + stimulation:leg + direction:leg + subject	0.0000795
stimulation + direction + stimulation:direction + leg + stimulation:leg + direction:leg + stimulation:direction:leg + subject	0.0000229

The winning model is the one with all main effects, with a Bayes factor of 1.9. We can compute the evidence for a particular effect by comparing this winning model with the best-fitting model that does not contain the effect. We can compute the evidence for absence of a particular effect by comparing the winning model with the best-fitting model that does contain the effect.

• STIMULATION. This effect can be quantified by comparing the Bayes factors of the first and the 3rd model : 6.9. This constitues moderate evidence for the presence of a stimulation effect.

• DIRECTION. This effect can be quantified by comparing the first and the 2nd model : 5.8. This constitues anecdotal evidence for the presence of a direction effect.

• Interaction of STIMULATION and DIRECTION. This term is not in the winning model, and first shows up in the fourth model. The Bayes factor for the comparison of both models is 11.3. This constitues moderate evidence for the absence of an interaction between stimulation and direction.
• Interaction of STIMULATION, and LEG. Bayes factor: 70.8. This constitutes moderate evidence for the absence of an interaction between stimulation and leg.

• Three-way interaction. Bayes factor: 8.4607710^{4}. This constitutes strong evidence for the absence of an interaction between stimulation, leg and direction.

Test against the full model

Another option for quantifying evidence for a particular effect is to compare the full model to a model where that effect is omitted (whichModels = top"). The full model is stimulation + direction + stimulation:direction + leg + stimulation:leg + direction:leg + stimulation:direction:leg + subject.

bfKanaiFull = anovaBF(latency~stimulation*leg*direction+subject, data = data.frame(latencyMedianBaselineLateral), whichModels="top", whichRandom = "subject", progress = FALSE, iterations = 100000) # compute Bayes Factors
bfKanaiFull

Bayes factor top-down analysis
--------------
When effect is omitted from stimulation + direction + leg + stimulation:direction + stimulation:leg + direction:leg + stimulation:direction:leg + subject , BF is...
[1] Omit direction:leg:stimulation : 3.531213  <U+00B1>2.8%
[2] Omit direction:leg             : 9.958357  <U+00B1>2.39%
[3] Omit leg:stimulation           : 12.17389  <U+00B1>2.93%
[4] Omit direction:stimulation     : 5.00078   <U+00B1>2.83%
[5] Omit leg                       : 0.5059081 <U+00B1>2.83%
[6] Omit direction                 : 7.274857  <U+00B1>2.61%
[7] Omit stimulation               : 5.96352   <U+00B1>3.1%

Against denominator:
  latency ~ stimulation + direction + leg + stimulation:direction + stimulation:leg + direction:leg + stimulation:direction:leg +     subject 
---
Bayes factor type: BFlinearModel, JZS

Removing the DIRECTION effect from the model yields a lower Bayes Factor, so including this effect improved the model. This thus constitues (some) evidence for the alternative: 1 \ 7.275 = 0.1. The evidence for a stimulation effect is also about the same: 0.2.

On the contrary, removing the interactions actually improves the model, so there is moderate evidence for the absence of interaction effects.

All in all, for most effects, the result is qualitatively similar to the approach of testing against the null-model. Only the evidence for absence of the three-way interaction is much, much smaller here.

Bayesian linear mixed effects matching Kanai - center saccades

bfKanaiCenter = anovaBF(latency~stimulation*leg*direction+subject, data = data.frame(latencyMedianBaselineCenter), whichModels="top", whichRandom = "subject", progress = FALSE, iterations = 100000) # compute Bayes Factors
bfKanaiCenter

Bayes factor top-down analysis
--------------
When effect is omitted from stimulation + direction + leg + stimulation:direction + stimulation:leg + direction:leg + stimulation:direction:leg + subject , BF is...
[1] Omit direction:leg:stimulation : 4.287449 <U+00B1>3.81%
[2] Omit direction:leg             : 7.782446 <U+00B1>2.85%
[3] Omit leg:stimulation           : 12.15958 <U+00B1>3.32%
[4] Omit direction:stimulation     : 5.243033 <U+00B1>4.65%
[5] Omit leg                       : 18.55041 <U+00B1>33.29%
[6] Omit direction                 : 1.008738 <U+00B1>46.93%
[7] Omit stimulation               : 6.57132  <U+00B1>2.11%

Against denominator:
  latency ~ stimulation + direction + leg + stimulation:direction + stimulation:leg + direction:leg + stimulation:direction:leg +     subject 
---
Bayes factor type: BFlinearModel, JZS

Bayesian linear mixed effects with baseline - saccade latency

Repeat the analysis when the baseline is not subtracted, and also for center saccades.

Lateral saccades

# keep only lateral saccades
latencyMedianLateral <- latencyMedianLegExcl %>%
  filter(type == "lateral") %>%
  select(-type)
  
bfLateral = anovaBF(latency~stimulation*leg*direction+subject, data = data.frame(latencyMedianLateral), whichModels="top", whichRandom = "subject", progress = FALSE, iterations = 100000) # compute Bayes Factors
bfLateral

Bayes factor top-down analysis
--------------
When effect is omitted from stimulation + direction + leg + stimulation:direction + stimulation:leg + direction:leg + stimulation:direction:leg + subject , BF is...
[1] Omit direction:leg:stimulation : 11.40207  <U+00B1>2.12%
[2] Omit direction:leg             : 28.10296  <U+00B1>1.9%
[3] Omit leg:stimulation           : 28.90428  <U+00B1>1.76%
[4] Omit direction:stimulation     : 3.637222  <U+00B1>1.88%
[5] Omit leg                       : 17.51425  <U+00B1>3.04%
[6] Omit direction                 : 2.407894  <U+00B1>1.74%
[7] Omit stimulation               : 0.5932669 <U+00B1>3.19%

Against denominator:
  latency ~ stimulation + direction + leg + stimulation:direction + stimulation:leg + direction:leg + stimulation:direction:leg +     subject 
---
Bayes factor type: BFlinearModel, JZS

When including the baseline, the evidence for null effects are stronger all around.

Center saccades

# keep only center saccades
latencyMedianCenter <- latencyMedianLegExcl %>%
  filter(type == "center") %>%
  select(-type)
  
bfCenter = anovaBF(latency~stimulation*leg*direction+subject, data = data.frame(latencyMedianCenter), whichModels="top", whichRandom = "subject", progress = FALSE, iterations = 100000) # compute Bayes Factors
bfCenter

Bayes factor top-down analysis
--------------
When effect is omitted from stimulation + direction + leg + stimulation:direction + stimulation:leg + direction:leg + stimulation:direction:leg + subject , BF is...
[1] Omit direction:leg:stimulation : 11.78923 <U+00B1>4.91%
[2] Omit direction:leg             : 19.688   <U+00B1>5.21%
[3] Omit leg:stimulation           : 31.49381 <U+00B1>7.63%
[4] Omit direction:stimulation     : 5.853747 <U+00B1>2.62%
[5] Omit leg                       : 44.50356 <U+00B1>2.83%
[6] Omit direction                 : 7.08912  <U+00B1>1.9%
[7] Omit stimulation               : 6.663299 <U+00B1>1.74%

Against denominator:
  latency ~ stimulation + direction + leg + stimulation:direction + stimulation:leg + direction:leg + stimulation:direction:leg +     subject 
---
Bayes factor type: BFlinearModel, JZS

Now the model improves all the time, no matter which effect is omitted.

Reaction time distributions

Split data in equal blocks

We want the same number of trials to go into each distribution, so we’ll split the 30-minute long post leg into two. Then we have 4 15-minute long data segments: pre, tDCS, post.1 and post.2.

# Group single-trial data in 4 legs instead of 3
groupDataBlocked <- groupData %>%
  mutate(leg = as.character(leg)) %>% #un-factor the leg column so we can change it
  mutate(leg=replace(leg, leg == "post", "post.1"), # split the "post" level into two
         leg=replace(leg, block > 3, "post.2")) %>%
  mutate(leg = factor(leg, levels = c("pre","tDCS","post.1","post.2"))) #re-factor

Fit the model

Next we fit an ex-gaussian distribution using the retimes package. One model is fit for each combination of:

• SUBJECT (S01, S02, etc.)
• STIMULATION (anodal, cathodal)
• LEG (pre, tDCS, post.1, post.2)
• TYPE (center, lateral)
• DIRECTION (left, right)

The ex-gaussian distribution is the convolution of a normal and an exponential distribution. It can be characterized with 3 parameters:

• \(\mu\) (mu): the mean of the normal distribution
• \(\sigma\) (sigma): the variance of the normal distribution
• \(\tau\) (tau): the rate parameter of the exponential distribution

We’ll do minimal pre-processing: keep all subjects, only discard negative RTs

groupDataFit <- groupDataBlocked %>%
  filter(latency > 0) %>% # throw out missing and negative RTs
  group_by(subject,stimulation,leg,type,direction) %>% # for every condition
  summarise(result = list(data.frame(t(attr(timefit(latency),"par"))))) %>% # [see below]
  ungroup() %>% # remove grouping information (otherwise unnest will not work)
  unnest() # unpack the list column, so each parameter gets its own column
  
# 1. fit the ex-gaussian to the latency data with the timefit() function
# 2. extract the fitted parameters from the resulting object with "attr()"
# 3. transpose ("t()") so every parameter is in a different column
# 4. convert to a data frame with one column per parameter
# 5. pack in a list so we have something of size 1 we can assign to the dataframe
# "result" is now a list-column, where each element is a dataframe with estimates of the 3 parameters
# Alternatively, create a function to fit the data, and then call it through dplyr::do
#myFunc <- function(x){
#  data.frame(t(attr(timefit(x),"par")))
#}
#groupDataFit <- groupDataBlocked %>%
#  group_by(subject,stimulation,leg,type,direction) %>% # for every condition
#  do((myFunc(.$latency)))

Inspect fitted parameters

ggplot(groupDataFit, aes(interaction(leg,stimulation), mu, color = subject)) +
  facet_grid(type~direction) + 
  stat_summary(fun.y = mean, fun.ymin = function(x) mean(x) - 2*sd(x), fun.ymax = function(x) mean(x) + 2*sd(x), aes(group = 1)) + # plot mean and 2*SD to watch for outliers; we need group = 1 because the x-axis is a combination of factors
  geom_jitter(width = 0.6) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

For center saccades some estimates of \(\mu\) are simply too fast (e.g. S09); they reflect premature eye movements. The distribution is also not a good fit in these cases.

For center saccades, it looks reasonable. S01 is more than 2 SD away from the mean, but the fit is excellent.

ggplot(groupDataFit, aes(interaction(leg,stimulation), sigma, color = subject)) +
  facet_grid(type~direction) +
  stat_summary(fun.y = mean, fun.ymin = function(x) mean(x) - 2*sd(x), fun.ymax = function(x) mean(x) + 2*sd(x), aes(group = 1)) +
  geom_jitter(width = 0.6) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

Center saccades appear to have larger \(\sigma\) than lateral. But if you inspect the distributions, the peak is actually much sharper, which should indicate lower variance. The problem is that the shape of many distributions does not conform well: many saccades are early and have little variance, but then there’s also a bump in the right-tail which causes the sigma-parameter to increase.

Values for the lateral distribution seem reasonable. 2 estimates of \(\sigma\) are too close to zero, others are too high, but the spread is not unreasonable. Inspecting these distributions also reveals presence of a “shoulder” right of the peak, which can either spread the distribution out or narrow it too much.

ggplot(groupDataFit, aes(interaction(leg,stimulation), tau, color = subject)) +
  facet_grid(type~direction) +
  stat_summary(fun.y = mean, fun.ymin = function(x) mean(x) - 2*sd(x), fun.ymax = function(x) mean(x) + 2*sd(x), aes(group = 1)) +
  geom_jitter(width = 0.6) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

Particularly for center saccades, \(\tau\) is many times estimated at 0. This means that the distribution is actually just a normal distribution. For center saccades that makes sense, as \(\tau\) is responsible for the longer right tail that characterizes the ex-gaussian distribution. When saccades are all super fast, this tail is not present.

Mu

Plot

muLinePlot <- ggplot(groupDataFit, aes(leg, mu, color = stimulation, shape = stimulation)) +         
  facet_grid(type ~ direction) +
  stat_summary(fun.y = mean, geom = "point", size = 3) +
  stat_summary(fun.y = mean, geom = "line", aes(group = stimulation), size = 1) +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", width = 0.3)
muLinePlot

Statistics

Lateral

muLateral <- groupDataFit %>%
  filter(type == "lateral") %>%
  select(-type)

modelMuLateral <- ezANOVA(data = data.frame(muLateral), dv = .(mu), wid = .(subject), within = .(stimulation, leg, direction), type = 3)
kable(modelMuLateral$ANOVA)

	Effect	DFn	DFd	F	p	p<.05	ges
2	stimulation	1	19	0.4962574	0.4896947		0.0018185
3	leg	3	57	4.0770965	0.0108122	*	0.0068828
4	direction	1	19	0.0862917	0.7721303		0.0007264
5	stimulation:leg	3	57	0.2632381	0.8515909		0.0003878
6	stimulation:direction	1	19	2.4797673	0.1318238		0.0016080
7	leg:direction	3	57	3.7885316	0.0150864	*	0.0027886
8	stimulation:leg:direction	3	57	1.1099797	0.3526354		0.0010857

bfMuLateral = anovaBF(mu~stimulation*leg*direction+subject, data = data.frame(muLateral), whichModels="top", whichRandom = "subject", progress = FALSE, iterations = 100000) # compute Bayes Factors
bfMuLateral

Bayes factor top-down analysis
--------------
When effect is omitted from stimulation + leg + direction + stimulation:leg + stimulation:direction + leg:direction + stimulation:leg:direction + subject , BF is...
[1] Omit direction:leg:stimulation : 10.44442 <U+00B1>5.83%
[2] Omit direction:leg             : 11.20198 <U+00B1>5.71%
[3] Omit direction:stimulation     : 2.92745  <U+00B1>5.11%
[4] Omit leg:stimulation           : 24.46828 <U+00B1>4.89%
[5] Omit direction                 : 5.742728 <U+00B1>5.39%
[6] Omit leg                       : 4.756499 <U+00B1>5.04%
[7] Omit stimulation               : 3.504989 <U+00B1>4.89%

Against denominator:
  mu ~ stimulation + leg + direction + stimulation:leg + stimulation:direction + leg:direction + stimulation:leg:direction + subject 
---
Bayes factor type: BFlinearModel, JZS

Center

muCenter <- groupDataFit %>%
  filter(type == "center") %>%
  select(-type)

modelMuCenter <- ezANOVA(data = data.frame(muCenter), dv = .(mu), wid = .(subject), within = .(stimulation, leg, direction), type = 3)
kable(modelMuCenter$ANOVA)

	Effect	DFn	DFd	F	p	p<.05	ges
2	stimulation	1	19	0.0001866	0.9892442		0.0000011
3	leg	3	57	1.0722467	0.3681380		0.0060317
4	direction	1	19	0.0003526	0.9852143		0.0000024
5	stimulation:leg	3	57	1.2683962	0.2938734		0.0043846
6	stimulation:direction	1	19	0.6514095	0.4295975		0.0010113
7	leg:direction	3	57	0.5453032	0.6533091		0.0020456
8	stimulation:leg:direction	3	57	0.6311947	0.5978988		0.0014300

bfMuCenter = anovaBF(mu~stimulation*leg*direction+subject, data = data.frame(muCenter), whichModels="top", whichRandom = "subject", progress = FALSE, iterations = 100000) # compute Bayes Factors
bfMuCenter

Bayes factor top-down analysis
--------------
When effect is omitted from stimulation + leg + direction + stimulation:leg + stimulation:direction + leg:direction + stimulation:leg:direction + subject , BF is...
[1] Omit direction:leg:stimulation : 11.44415 <U+00B1>3.15%
[2] Omit direction:leg             : 19.57465 <U+00B1>2.23%
[3] Omit direction:stimulation     : 4.617263 <U+00B1>2.46%
[4] Omit leg:stimulation           : 12.14351 <U+00B1>2.29%
[5] Omit direction                 : 8.223817 <U+00B1>3.08%
[6] Omit leg                       : 17.60705 <U+00B1>2.35%
[7] Omit stimulation               : 8.096295 <U+00B1>2.49%

Against denominator:
  mu ~ stimulation + leg + direction + stimulation:leg + stimulation:direction + leg:direction + stimulation:leg:direction + subject 
---
Bayes factor type: BFlinearModel, JZS

Sigma

Plot

sigmaLinePlot <- ggplot(groupDataFit, aes(leg, sigma, color = stimulation, shape = stimulation)) +         
  facet_grid(type ~ direction) +
  stat_summary(fun.y = mean, geom = "point", size = 3) +
  stat_summary(fun.y = mean, geom = "line", aes(group = stimulation), size = 1) +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", width = 0.3)
sigmaLinePlot

Statistics

Lateral

sigmaLateral <- groupDataFit %>%
  filter(type == "lateral") %>%
  select(-type)

modelSigmaLateral <- ezANOVA(data = data.frame(sigmaLateral), dv = .(sigma), wid = .(subject), within = .(stimulation, leg, direction), type = 3)
kable(modelSigmaLateral$ANOVA)

	Effect	DFn	DFd	F	p	p<.05	ges
2	stimulation	1	19	0.3070893	0.5859358		0.0020213
3	leg	3	57	4.8020840	0.0047376	*	0.0205278
4	direction	1	19	0.4964608	0.4896068		0.0019499
5	stimulation:leg	3	57	0.3144879	0.8148281		0.0015190
6	stimulation:direction	1	19	0.8417414	0.3703991		0.0012920
7	leg:direction	3	57	1.4989067	0.2245956		0.0037770
8	stimulation:leg:direction	3	57	1.5867054	0.2025735		0.0049500

bfSigmaLateral = anovaBF(sigma~stimulation*leg*direction+subject, data = data.frame(sigmaLateral), whichModels="top", whichRandom = "subject", progress = FALSE, iterations = 100000) # compute Bayes Factors
bfSigmaLateral

Bayes factor top-down analysis
--------------
When effect is omitted from stimulation + leg + direction + stimulation:leg + stimulation:direction + leg:direction + stimulation:leg:direction + subject , BF is...
[1] Omit direction:leg:stimulation : 5.269379  <U+00B1>2.38%
[2] Omit direction:leg             : 13.71165  <U+00B1>5.71%
[3] Omit direction:stimulation     : 4.112387  <U+00B1>2.31%
[4] Omit leg:stimulation           : 24.85324  <U+00B1>11%
[5] Omit direction                 : 4.874215  <U+00B1>2.83%
[6] Omit leg                       : 0.4419473 <U+00B1>10.42%
[7] Omit stimulation               : 4.648076  <U+00B1>3.3%

Against denominator:
  sigma ~ stimulation + leg + direction + stimulation:leg + stimulation:direction + leg:direction + stimulation:leg:direction + subject 
---
Bayes factor type: BFlinearModel, JZS

Center

sigmaCenter <- groupDataFit %>%
  filter(type == "center") %>%
  select(-type)

modelSigmaCenter <- ezANOVA(data = data.frame(sigmaCenter), dv = .(sigma), wid = .(subject), within = .(stimulation, leg, direction), type = 3)
kable(modelSigmaCenter$ANOVA)

	Effect	DFn	DFd	F	p	p<.05	ges
2	stimulation	1	19	0.6283303	0.4377549		0.0037630
3	leg	3	57	5.2100287	0.0030012	*	0.0195091
4	direction	1	19	3.1389212	0.0924826		0.0211137
5	stimulation:leg	3	57	0.5290551	0.6641446		0.0026847
6	stimulation:direction	1	19	0.0135394	0.9085893		0.0000365
7	leg:direction	3	57	0.6716000	0.5729879		0.0032568
8	stimulation:leg:direction	3	57	0.7492325	0.5273056		0.0019792

bfSigmaCenter = anovaBF(sigma~stimulation*leg*direction+subject, data = data.frame(sigmaCenter), whichModels="top", whichRandom = "subject", progress = FALSE, iterations = 100000) # compute Bayes Factors
bfSigmaCenter

Bayes factor top-down analysis
--------------
When effect is omitted from stimulation + leg + direction + stimulation:leg + stimulation:direction + leg:direction + stimulation:leg:direction + subject , BF is...
[1] Omit direction:leg:stimulation : 10.06561   <U+00B1>3.98%
[2] Omit direction:leg             : 15.74758   <U+00B1>3.55%
[3] Omit direction:stimulation     : 5.793119   <U+00B1>3.68%
[4] Omit leg:stimulation           : 17.38659   <U+00B1>3.72%
[5] Omit direction                 : 0.05443641 <U+00B1>3.56%
[6] Omit leg                       : 1.08491    <U+00B1>3.98%
[7] Omit stimulation               : 3.826635   <U+00B1>10.66%

Against denominator:
  sigma ~ stimulation + leg + direction + stimulation:leg + stimulation:direction + leg:direction + stimulation:leg:direction + subject 
---
Bayes factor type: BFlinearModel, JZS

Tau

Plot

tauLinePlot <- ggplot(groupDataFit, aes(leg, tau, color = stimulation, shape = stimulation)) +         
  facet_grid(type ~ direction) +
  stat_summary(fun.y = mean, geom = "point", size = 3) +
  stat_summary(fun.y = mean, geom = "line", aes(group = stimulation), size = 1) +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", width = 0.3)
tauLinePlot

Statistics

Lateral

tauLateral <- groupDataFit %>%
  filter(type == "lateral") %>%
  select(-type)

modelTauLateral <- ezANOVA(data = data.frame(tauLateral), dv = .(tau), wid = .(subject), within = .(stimulation, leg, direction), type = 3)
kable(modelTauLateral$ANOVA)

	Effect	DFn	DFd	F	p	p<.05	ges
2	stimulation	1	19	2.4821341	0.1316505		0.0108311
3	leg	3	57	2.7975288	0.0481938	*	0.0178416
4	direction	1	19	0.0020449	0.9644036		0.0000099
5	stimulation:leg	3	57	0.5265115	0.6658504		0.0018047
6	stimulation:direction	1	19	0.0180463	0.8945494		0.0000357
7	leg:direction	3	57	2.0908871	0.1114990		0.0045650
8	stimulation:leg:direction	3	57	0.3432422	0.7941361		0.0011038

bfTauLateral = anovaBF(tau~stimulation*leg*direction+subject, data = data.frame(tauLateral), whichModels="top", whichRandom = "subject", progress = FALSE, iterations = 100000) # compute Bayes Factors
bfTauLateral

Bayes factor top-down analysis
--------------
When effect is omitted from stimulation + leg + direction + stimulation:leg + stimulation:direction + leg:direction + stimulation:leg:direction + subject , BF is...
[1] Omit direction:leg:stimulation : 11.00158  <U+00B1>5.51%
[2] Omit direction:leg             : 10.05369  <U+00B1>5.53%
[3] Omit direction:stimulation     : 5.436445  <U+00B1>5.41%
[4] Omit leg:stimulation           : 31.78977  <U+00B1>37.27%
[5] Omit direction                 : 7.690549  <U+00B1>5.6%
[6] Omit leg                       : 0.6883419 <U+00B1>5.7%
[7] Omit stimulation               : 0.3571577 <U+00B1>5.46%

Against denominator:
  tau ~ stimulation + leg + direction + stimulation:leg + stimulation:direction + leg:direction + stimulation:leg:direction + subject 
---
Bayes factor type: BFlinearModel, JZS

Center

tauCenter <- groupDataFit %>%
  filter(type == "center") %>%
  select(-type)

modelTauCenter <- ezANOVA(data = data.frame(tauCenter), dv = .(tau), wid = .(subject), within = .(stimulation, leg, direction), type = 3)
kable(modelTauCenter$ANOVA)

	Effect	DFn	DFd	F	p	p<.05	ges
2	stimulation	1	19	0.9421655	0.3439170		0.0049757
3	leg	3	57	0.9034300	0.4452132		0.0068505
4	direction	1	19	0.4129620	0.5281491		0.0016242
5	stimulation:leg	3	57	2.3352550	0.0833949		0.0095634
6	stimulation:direction	1	19	0.6623347	0.4258193		0.0014819
7	leg:direction	3	57	1.1465225	0.3381913		0.0067472
8	stimulation:leg:direction	3	57	0.1164547	0.9501026		0.0004604

bfTauCenter = anovaBF(tau~stimulation*leg*direction+subject, data = data.frame(tauCenter), whichModels="top", whichRandom = "subject", progress = FALSE, iterations = 100000) # compute Bayes Factors
bfTauCenter

Bayes factor top-down analysis
--------------
When effect is omitted from stimulation + leg + direction + stimulation:leg + stimulation:direction + leg:direction + stimulation:leg:direction + subject , BF is...
[1] Omit direction:leg:stimulation : 13.92171 <U+00B1>3.97%
[2] Omit direction:leg             : 8.539389 <U+00B1>3.44%
[3] Omit direction:stimulation     : 4.293162 <U+00B1>3.61%
[4] Omit leg:stimulation           : 5.506957 <U+00B1>5.76%
[5] Omit direction                 : 5.55446  <U+00B1>3.29%
[6] Omit leg                       : 18.09544 <U+00B1>3.77%
[7] Omit stimulation               : 2.643075 <U+00B1>3.74%

Against denominator:
  tau ~ stimulation + leg + direction + stimulation:leg + stimulation:direction + leg:direction + stimulation:leg:direction + subject 
---
Bayes factor type: BFlinearModel, JZS